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1 Introduction 

Variational integrators (Vis) [10, 7, 8^ are a general class of time integration methods for Hamiltonian systems 
whose construction guarantees certain highly desirable properties. Instead of directly discretizing the smooth 
equations of motion of a system, the variational approach asks that we instead step back and discretize 
the system's Lagrangian. By analogy to Hamilton's Least Action Principle, we may then form a discrete 
action and seek paths which extremize it, yielding discrete Euler-Lagrange equations from which discrete 
equations of motion are readily recovered. As a consequence of this special, more principled construction, 
variational integrators are guaranteed to satisfy a discrete formulation of Noether's Theorem [T2], and as a 
special case conserve linear and angular momentum. Vis are automatically symplectic [2] ; while they do not 
necessarily conserve energy, conservation of the symplectic form assures no-drift conservation of energy over 
exponentially many time steps [2|. 

Mechanical systems are almost never uniformly stiff. Different potentials have different stable time step 
requirements, and even for identical potentials this requirement depends on element size, since finer elements 
can support higher-energy modes than coarser elements. Any global time integration scheme cannot take 
advantage of this variability, and instead must integrate the entire system at the globally stiffest time step. 
Suppose the system can be triangulated into elements such that each force acts entirely within one element. 
Then asynchronous variational integrators [6] generalize Vis by allowing each element to have its own, 
independent time step. Coarser elements can then be assigned a slower "clock," and finer elements a faster 
one, so that relatively few very fine elements do not as significantly degrade the overall performance of 
integrating the system. AVIs retain all of the properties of variational integrators mentioned above, except 
for symplecticity. However, AVIs instead preserve an analogous multisymplectic form, and it has been shown 
experimentally that preservation of this form likely induces the same long-time good energy behavior that 
characterize symplectic integrators [6]. 

The published proof of multisymplecticity assumes that the potentials are of an "elastic type," i.e., spec- 
ified by volume integration over the material domain, an assumption violated by interaction-type potentials. 
We extend the proof, showing that AVIs remain multisymplectic under relaxed assumptions on the type of 
potential (S}2j). The modified proof allows for interaction potentials of the kind needed for contact mechan- 
ics (i.e., penalty forces). The extended theory thus enables the simulation of mechanical contact in elastica 
(such as thin shells) and multibody systems (such as granular materials) with no drift of conserved quantities 
(energy, momentum) over long run times, using the algorithms in [3J. 

We conclude with data from numerical experiments measuring the long time energy behavior of simulated 
contact, comparing the method built on multisymplectic integration of interaction potentials to recently 
proposed methods for thin shell contact 

2 Variational Integrators 

We begin with a background on variational integration and symplectic structure [21 [8l I12j . 
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Let 7(t) be a piecewise-regular trajectory through configuration space Q, and j(t) = ^l{t) be the 
configurational velocity at time t. For simplicity we shall assume that the kinetic energy of the system T 
depends only on configurational velocity, and that the potential energy V depends only on configurational 
position, so that we may write the Lagrangian L at time t as 

L(q,q)=T(q)-V(q). (1) 

Then given the configuration of the system qg at time to and qf at tf, Hamilton's principle [1] states 
that the trajectory of the system j(t) joining j(to) — qo and 7(t/) = qf is a stationary point of the action 
functional 

S(l)= f 1 L[y(t),j(t)]dt 



with respect to taking variations £7 of 7 which leave 7 fixed at the endpoints to, tf. In other words, 7 
satisfies 

dS(j) ■ 67 = 0. (2) 
Integrating by parts, and using that Sj vanishes at to and ti, we compute 

rtf ( dL , , „ dL, , „.\ . f tf f dV , , d 2 T , 



dS(j) ■ Sj 



tt \"i / Jt 

Since this equality must hold for all variations $7 that fix 7's endpoints, we must have 

dV d 2 T 

^(7) + 7 ^(^ = °' ( 3 ) 

the Euler- Lagrange equation of the system. This equation is a second-order ordinary differential equation, 
and so has a unique solution 7 given two initial values 7(^0) fmd 7(^0) ■ 

2.1 Symplecticity 

The flow 6 S : ['y{t), A f(t)] [j(t + s),j(t + s)] given by §3§ has many structure-preserving properties; in 
particular it is momentum-preserving, energy-preserving, and symplectic [5]. To see this last property, for 
the remainder of this section we restrict the space of trajectories to those that satisfy the Euler-Lagrange 
equations. For such trajectories, and relaxing the requirement that ^7 fix the endpoints of 7, we have 



dT 

dS{i) -6-f = — [iTq(q, q)\ ■ Sj 



(4) 



where itq is projection onto the second factor. 

Since initial conditions (q, q) are in bijection with trajectories satisfying the Euler-Lagrange equation, 
such trajectories 7 can be uniquely parameterized by initial conditions [7(to), 7(^0)]- F° r the remainder of 
this section we also restrict variations S"f to first variations: those variations in whose direction 7 continues 
to satisfy the Euler-Lagrange equations. These are also parameterized by variations of the initial conditions, 
(5q,5q). For conciseness of notation, we will write v(t) = (j(t),j(t)) and 5v(t) — [Sj(t), Sjit)]; using this 
notation we write the above two facts as v(t) = Qt-t v(to) and 8v(t) = Q t -t 0tt o~v(to)- The action ([T]), a 
functional on trajectories 7, can also be rewritten as a function Si of the initial conditions, 

Si{qA)= I' ° L[Q t (q,q)}dt, 
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so that 

ds( 7 ) . 5 7 = d$ [i/(t )] • <yi/(t ) 

Substituting all of these expressions into Q , we get 

'ar 



dSi [u(t )] ■ Sv(to) 



tf 



dT 



o 7r 4 [8 t _ io i/(i )] d<? • <M*) 



'<9T 

= ( gT- ° TTg ) [8 t -t I/(*o)] d<? • ©t-t *M*o) 



where #l is the one-form o 7r<j^ dg. Since dSj is exact, 



d 2 Si = = 6 



t/-to 



'de L -de L , 



so since to and tf are arbitrary, <3* s d0L = ddL for arbitrary times s, and preserves the so-called symplectic 
form ddh- 

2.2 Discretization 

Discrete mechanics [TTJ [101 [SI HI H] describes a discretization of Hamilton's principle, yielding a numerical 
integrator that shares many of the structure-preserving properties of the continuous flow S . Consider 
a discretization of the trajectory 7 : [to,t/] — > Q by a piecewise linear trajectory interpolating n points 
q = {<7o:9i 5 ■ ■ with go = 7(^0) and q n -i = "/{tf), where the discrete velocity g,-|_i/2 on the segment 

between and g^+i is 



<?i+l/2 



ft = 



ft ' n — 1 

We seek an analogue of ([3]) in this discrete setting. To that end, we formulate a discrete Lagrangian 

' qb - q a 



Ld{q a ,qb) = T — - V(q b 



(5) 



and discrete action 



Sd(q) = hL d (q u q i+ i). (6) 

i=0 

Motivated by ([2]), we impose a discrete Hamilton's principle: 

dSM) • 5q = 

for all variations <5q = {<5qoj ■ ■ ■ j &?n— 1} that fix q at its endpoints, i.e. , with Sq^ — 5q n _i = 0. For 
ease of notation, we define versions of the kinetic and potential energy terms in ([5]) that depend on (q a ,q b ) 
instead of (q, q): 



Td{q a ,qb) = T 



qb - q a 



V d (q a ,qb) = V(q b ) 



rpi, \ dT f q b -q a 
dV 

Vd(q a ,qb) = -^-(qb)- 
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Then 



n-2 



dS d (q) ■ Sq = 2J h (DiL d (q il q i+ i) ■ Sq { + D 2 L d (qi,q i+1 ) ■ Sq i+ i) 

dq 



i=0 
n-2 



]^ h (~T T d(H>1i+l) ■ 5q % + ^-T d (gi, q i+1 ) ■ Sq i+1 - ^-{q l+ i) ■ %- 



i=0 



= 7d(g n _2,g n _i) • (5q„-i - T' d {q ,qi) ■ Sq - h—{q n ^x) ■ Sq n _i 
+ Y1 ( T d(qi-i,qi) -T' d {q u q i+ i) J • % 

= X] (^ T dfe-i,9i) - ^d(%,g i+ i) - h— (q t )J ■8q i = 0. 

Since is unconstrained for 1 < i < n — 2, we must have 

dT dT dV 

- - -/i-T^-(gi), i = l,...,n-2, (7) 

the discrete Euler-Langrange equations of the system. 

Unlike in the continuous settings, the discrete Euler-Lagrange equations do not always have a unique 
solution given initial values qo and q±. We therefore assume in all that follows that T d and V d are of a form 
so that ([7]) gives a unique q^+i given qi and (7^-1 — this assumption always holds, for instance, in the typical 
case where T d is quadratic in q. Then the discrete Euler-Lagrange equations give a well-defined discrete flow 

F : {q l -i,q % ) i-> (fc+i), 

which recovers the entire trajectory from initial conditions, in perfect analogy to the continuous setting. 
2.3 Symplecticity of the Discrete Flow 

We now would like a symplectic form preserved by F, just as d0L is preserved by 0. As in the continuous 
setting, we restrict trajectories q to those that satisfy the discrete Euler-Lagrange equations, and restrict 
variations to first variations (and relax the condition that these variations vanish at the endpoints) , yielding 

dV 

dS d (q) ■ Sq = T d (g n _ 2 , q n -x) • 5q n -i - T d (q ,qi) ■ 5q - h—(q n -i) ■ 5q n -i- 

We denote by F k the discrete flow F composed with itself k times, or k "steps" of F. We remark again 
that all q satisfying ([7]) can be parameterized by initial conditions Vq = (qa,qi), and first variations by 
5vq = (Sqo, Sqi), so that we can rewrite the discrete action as 

n-2 



Sid(vo) = ^2hL d (F i Vti). 



i=0 
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Putting together all of the pieces, 



dS id (v ) ■ 5v = dS d (q) ■ Sq 

dV 

= T d(Qn-2,q n -i) ■ Sq n -i - T d (q , qi) ■ 5q - h—(q n -i) ■ 6q n ~i 




'a.=q-n.-2, q b =q n -i 



- T d(la, qb)dq a ■ (Sqo,Sqi) 

9a =90, 96=91 

= [T' d (F n - 2 v Q ) - hV'(F n ~ 2 i> )] dq b ■ F n -\5v - T' d { Vo )dq a ■ 8v 
= °U-* Vo ■ F n - 2 Jv + 9- ■ 5v Q 
= (F n - 2 *9+) ■ 5u + 9- ■ 5v . 

for the indicated two-forms 9 + and 9~ . Since d(hL d ) = 9 + + 9~ , d 2 (hL d ) — = d9 + + d9~ . Moreover the 
intial conditions vq are arbitrary, hence 

d 2 S ld = = F n - 2 *d9 + + d0- = -F n - 2 *d0- + d6~, 



so 

dS~ =F n - 2 *d0-. 



Since n is arbitrary, we conclude that the discrete flow F preserves the symplectic form d9~ . Using backwards 
error analysis, it can be shown that this geometric property guarantees that integrating with F introduces 
no energy drift for a number of steps exponential in h [2], a highly desirable property when simulating 
molecular dynamic or other Hamiltonian systems whose qualitative behavior is substantially affected by 
errors in energy 



3 Asynchronous Variational Integrators 

In section 12.21 we formulated an action functional ((6]) as the integration of a single discrete Lagrangian 
over a single time step size h. Such a construction is cumbersome when modeling multiple potentials of 
varying stiffnesses acting on different parts of the system: to prevent instability we are forced to integrate 
the entire system at the resolution of the stiffest force. Given a spatial triangulation T = {K} of the 
system, asynchronous variational integrators (AVIs) , introduced by Lew et al. [6] , are a family of numerical 
integrators, derived from a discrete Hamilton's principle, that support integrating potentials on different 
triangles at different time steps. In the exposition that follows, we follow the arguments set forth by Lew 
et al., but depart at times from the notation used in their work. Although the additional notation and 
indices introduced herein are initially cumbersome, they will allow for a relatively easy transition to the 
triangulation- free setting in Section |U 

Instead of a global discrete Lagrangian, we instead imbue each triangle K with a local discrete Lagrangian 

L$(£, tf) = 1^ T K [q K (t)] dt - h K V K {q«\ 

where T K and V K are the elemental kinetic and potential energies on triangle K, respectively, h K = tb—t a is 
the elemental time step, and q K (t), the elemental velocity at time t, is left imprecise for the moment. We no 
longer assume that velocity is constant between times t a and if, — this would only be true if for every potential 
on a triangle adjacent to K, no multiple of its time step lies between t a and tf,, which is not necessarily the 
case — so unlike for the discrete Lagrangian (|5|), here we cannot explicitly integrate the kinetic energy term. 
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For this reason we now write the Lagrangian as an integrated quantity, instead of deferring the integration 
to inside the action. 

Each triangle is only concerned with certain moments in time — namely, integer multiples of h K — and 
these moments are inconsistent across triangles. We therefore subdivide time in a way compatible with all 
triangles: for a r-length interval of time, we define 

\r/h K \ 

= U U 

KeT j=o 

That is, 5(r) is the set of all integer multiples less than r of all elemental time steps. H can be ordered, 
and in particular we let be the {i + l)-st least clement of S. If n is the cardinality S, wc then discrctizc 
a trajectory of duration t by linearly interpolating intermediate configurations (/o, . . . , Qn— i> 

where q; is 

the configuration of the system at time We discretize velocity as qk+1/2 = g(fc+i)— g(fc) on tne se g men t 
of the trajectory between qk and ft+i. We now need to write a global action functional of these trajectories 
that sums the above elemental Lagrangians, which we do in the natural way: 

[T/h K \ 

5 A Vl(q)= E E L d (<?jW+l)- (8) 

As before, we consider variations <5q = {<5g , • • ■ , &?n-i} with 5g = <5<7n-i = 0, and impose Hamilton's 
principle, 

dSAvi(q) • 5q = 0. 

To avoid becoming bogged down in notation, we let co K (j) — £~ 1 (jh K ) — that is, u> maps local time indices 
for K to global indices into S — and will write qj interchangeably for nKQj, the restriction of the (global) 
configuration qj to an elemental configuration on K. Then 

wq) = E IE 1 L * (#'#+0 + f TK * 

(\jlh K \-\ 

= E E L * (^o>^o+d) + / r K [«*(*)] * 
= E E / rK [«*(*)] * - ^'(^w+d) + / T * [«*(*)] * 



Thus, writing 

nKflb-Qa\ rrK' r „ „_+ \ _ f Qb ~ Qa 

tb - t a 



T d K (q a ,q b ,t a ,t b ) = T K (f^) Tf'iq^ta,^) = ^~ ( 



( M ) = V^fe) V d K (q a ,q b ) = (<?&), 
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we have 



/n-2 i \ Vr/h K \ dyK 

dS A vi(q) -<5q= X! E T * fe' 9k+x, Z(k), £,(k + 1)] • (<5q fc+ i - Sq k ) - ^ h K ' 

kgT \k=o J KeT j=i q 



J2 (rf [qn-2,qn-i,t(n - 2),f(n - 1)] • «5<z„-i - if' [gb,9i,f(0),f(l)] ■ <fyo) 
r <eT 

n-2 

^ E E( T i r '[*- l '*^( fc - 1 )'^)]- T i r '[*'9 fc + l ^( fc )'^ + 1 )])-^ 



E E hK -g-jr(^ K u)) - S( i^ k u) 

KeT j=i q 



n-2 



E E ( T d'{Vk-l,qk,ak-l)A(k)}-Tf[q k ,q k+1 ,ak),ak + l)})-6q k 



k=l KeT 

™-2 p, v K 

"E E ^^r(ft)-*«. 
fc=ifcK|^(fe) 

where we abuse the notation /i^|m to mean, "all elemental time steps h which evenly divide m." Writing 
the total kinetic energy of the system ^2 Kel -T K as T to t, for AVIs we recover the discrete Euler-Lagrange 
equations 

(&-1/2) = - hK %K^)- (9) 

These equations are similar to those we derived for synchronous variational integrators ([7]), except that only 
a subset of potentials V\ contribute during each time step. As in the synchronous case, if, as is typical, T to t 
is quadratic in q, the system (jH]) gives rise to an explicit numerical integrator that is particularly easy to 
implement in practice. 

3.1 Multisymplecticity 

The right hand side of © depends on £(£;), and so unlike ([7]), the Euler-Lagrange equations for AVIs are time 
dependent, and do not give rise to a stationary update rule F(qi-i, qi) 1— » (gj,Qj + i). Instead, we consider 
the total, time-dependent flow F k (q ,qi) ^ (q k7 qk+i)- Once again, we parameterize trajectories satisfying 
© by vq = (<7o,<7i), and first variations by 8vq = (Sqo,Sqi). Restricting ourselves to such trajectories and 
variations, we rewrite the action ([5]) as 

(n-2 \r/h K \-l 
Y,m + l)-Z(k)}T?(F k (vo),ttk),!;(k + l)) - h K V d K (F- K ^- 1 M) 
k=0 j=0 
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Then 



dS iAV i(v) ■ 5v = dS AV i(q) ■ <5q 

= E ( T d' [<7n-2,<7„-i,£(n ~ 2), f(n - 1)] ■ <5<z„-i - T? [q , Ql , £(0), £(1)] ■ <%, 



E 



9g 



A' 



(<?n-l) • <5<?n- 



/i K |£(n-l) 

£ ( T f [^»-a( Wb ), € („-2),e(n-l)] •Jfe-x-Tf' [^o^(O), e(l)]-^o 
A'er 

- E W 

h K \H(n-l) 



= (0- + F n - 2 *0+U ■ 5v Q 
for one-forms 6~ and 6 + . Once again we have that 



= d^S = d6~ + F n ~ 2 *d9 + , 



(10) 



but unlike when our action was a sum of Lagrangians, from the multisymplectic form formula (|10[) we have 
no way of relating d0~ to d6 + , and thus do not recover symplectic structure preservation. Nevertheless, Lew 
et al. [6] conjecture that this multisymplectic structure leads to the good energy behavior observed for AVIs. 



4 Triangulat ion-Free AVIs 

The above formulation of AVIs assumed a spatial triangulation over which we defined distinct, local La- 
grangians. We now present a simple extension that supports potentials with arbitrary, possibly non-disjoint 
spatial stencil. 

Let {V 1 } be potentials with time steps h % . As in AVIs, for trajectories of duration r we define the set of 
times 

H(r) =[JU 3h\ 

V* 3=0 

the smallest set of times compatible with the time steps of all of the potentials. Again, let E have cardinality 
n, be the (i + l)-th least element of S, and = £ (j/i*)- Then, for T(q) the kinetic energy of the 

entire configuration, T d (q a , q b , t a , t b ) = T (f^5ff") > an d T' d (q a ,q b , t a , t b ) = ^? ( f^5ff ) j we write tne action 

n-2 Y T / h '\ 

sm) = E kcj + 1) - m T d faqj+uwuv + 1)] - E E ^n^). 

3=0 V 3=1 

We have made no attempt to define a Lagrangian pairing the kinetic and potential energy terms; we will see 
that an action defined this way still leads to a multisymplectic numeric integrator. 

To that end we impose dS g (q) • 5q = for variations with 5qo = 5q n -i = 0. Then we rewrite S g as 

n — 2 n— 1 

S g (c l ) = Y,lti3 + l)-a3)}T d lq j ,q j+ uZ(3),Z(3 + l)]-J2 E ^fo) 

3=0 3=lft*|e(3) 
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so that 



n-2 n-1 „ 

dS g (q) ■5q = J2^ WUV + 1)] • (<%+! - <%) - E E h V fe) ' ^ 



7^ [g„- 2 , <Z„-i,£(n - 2), £(n - 1)] • 5<? n -i - ^ [go, gi, e(0), C(l)] • 
i9g 



^IJCn-l) 
n-2 



3=1 \ h'lfW) 9 



n-2 



£ I ^ g,-, e(j - 1), CO')] - ^ fe, £0" + 1)] - E ] ' <%■ 

3=1 V 9 



The Euler-Lagrange equations are then 

dT 



M|?(fe) 



exactly the same as the Euler-Lagrange equations © for ordinary AVIs. Triangulation-free AVIs can thus 
be integrated in exactly the same manner as ordinary AVIs, for instance, by using the algorithm presented 
by Lew et al. [6]. 

4.1 Multisymplecticity 

To show that triangulation-free AVIs still satisfy the multisymplectic form formula (|10|). we follow the 
derivation for multisymplecticity of ordinary AVIs. Replacing X^ksT with Td in Section f3.ll an identical 
calculation shows triangulation-free AVIs satisfy (flQ]) . 



5 Sphere-plate Impact 

Our triangulation-free multisymplectic formulation supports interaction potentials of the kind needed for 
contact response, and we expect any such method to exhibit the good energy behavior associated with 
multisymplectic integrators. As a numerical experiment of this behavior, we simulated the impact of a 
spherical shell with a thin plate, as described in Cirak and West's article on Decomposition Contact Response 
(DCR) 1 , using the Asynchronous Contact Mechanics (ACM) framework [3] built on triangulation-free AVIs. 
A sphere of radius 0.125 approaches a plate of radius 0.35 with relative velocity 100. Both the sphere and 
the plate have thickness 0.0035. The time steps of our material forces (stretching and bending) are 10~ 7 
(the same as those chosen by Cirak and West.) 

Figure [T] compares energy over time when this simulation is run using both our ACM and DCR. Using 
ACM there is no noticeable long-term drift. Closely examining the energy data produced by ACM reveals the 
high-frequency, low-amplitude, qualitatively-negligible oscillations characteristic of symplectic integrators. 
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This work was supported in part by the NSF (MSPA Award No. IIS-05-28402, CSR Award No. CNS-06- 
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Figure 1: Total energy over time of a thin sphere colliding against a thin plate, simulated using asynchronous 
contact mechanics [3j (red) compared to data provided for decomposition contact response [I] (dark blue). 
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